#'
#' K-S method of niche breath and niche overlap estimation
#'  
#' @param sp1 the environmental variables of each individuals of sp1
#' @param sp2 the environmental variables of each individuals of sp2
#' @param nx number of bins used to estimated the emprical cumulative sample distribution.
#'
#' @details 
#' This is a continuous measurement of niche breadth and overlap. 
#' It defined as follows. Given two species A and B and the habitat variable,
#' e.g. altitudes, occupied by each individual, the altitudes for each species
#' were sorted into ascending order. It then calculated the cumulative sample
#' distribution P_A(x) and P_b(x), the proportion of stems of species A and species B
#' whose altitude is less than or equal to x. Niche overlap was caculated as twice the
#' area between the curve formed by the pairs of points (P_A,P_B) and either the 
#' P_A or P_B axis, whichever was smaller, resulting in a measure ranging between
#' 0 and 1. If two species did not overlap in altitudinal range, then the curve would 
#' lie along either the P_A axis and the line P_A=1 or along the P_B axis and the line P_B=1
#' and the niche overlap would be zero. If two species have identical altitudinal ranges,
#' P_A and P_B lie along the 45 degree line and the area between curve and either axis is
#' one-half which leads to a calculated niche overlap of one. 
#' 
#' Niche breadth for individual species was calcualted using a similar procedure. 
#' In this case, each species' distribution of altitudes was tested against a random sample
#' of 1000 points in the plot.
#' 
#' @references 
#' Potts, M. D., et al. (2004). "Habitat heterogeneity and niche structure of trees in two tropical rain forests." Oecologia 139: 446-453.
#'
#' @examples 
#' sp1=rnorm(100,mean=100)
#' sp2=rnorm(100,mean=0)
#' sp3=rnorm(100,mean=0)
#' 
#' nicheKS(sp1,sp2)
#' nicheKS(sp2,sp3)
#' 
#' @export
nicheKS=function(sp1,sp2=NULL,nx=10001){

  if(is.null(sp2))
    stop("can not calculate the niche overlap of a given species, the second species is missing.")
  
  x=seq(min(c(sp1,sp2)),max(c(sp1,sp2)),length.out=nx)
  p1=ecdf(sp1)(x)
  p2=ecdf(sp2)(x)
  stepx=1/(nx-1)
  nozero=p1!=0 & p2!=0
  re=sum(stepx*p2[nozero][-nx])
  if(re>0.5)
    re=1-re
  return(re*2)
}

